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Abstract 


We propose PALPS, a Process Algebra with Locations for Popula- 
tion Systems. PALPS allows us to produce spatially-explicit individual- 
based ecological models and to reason about their behavior. PALPS 
has two abstraction levels: At the first level, we may define the be- 
havior of an individual of a population and, at the second level, we 
may specify a system as the collection of individuals of various species 
located in space. In PALPS, the individuals move through their life 
cycle while changing their location and interact with each other in 
various ways such as predation, infection or mating. Furthermore, we 
propose a translation of a subset of PALPS into the probabilistic model 
checker PRISM. We illustrate our framework via models of dispersal in 
metapopulations and by applying PRISM on PALPS models for verifying 
temporal logic properties and conducting reachability and steady-state 
analysis. 

Keywords: Process calculi, ecology, probabilistic model checking, 
spatially-explicit individual-based models 


1 Introduction 


Population ecology is a subfield of ecology that deals with the dynamics of 
species populations and how these populations interact with the environ- 
ment. Its main aim is to understand how the population sizes of species 
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that live together in groups change over time and space. It has been of spe- 
cial interest to conservation scientists and practitioners who are interested 
in predicting how species will respond to specific management schemes and 
in guiding the selection of reservation sites and reintroduction efforts, e.g. 
[21, 36]. To make such predictions, scientists have been constructing models 
of ecosystems. These models are abstract representations of the systems in 
question which are subsequently studied to gain understanding of the real 
systems. Models enable researchers to simulate large-scale experiments that 
would be too costly or unethical to perform on a real ecosystem. They also 
enable the simulation of ecological processes over long periods of time. 


One of the main streams of today’s theoretical ecology is the individual- 
based approach to modeling population dynamics. In this approach, the 
modeling unit is that of a discrete individual and a system is considered as 
the composition of individuals and their environment. The concept of an 
individual does not always need to coincide with that of an individual or- 
ganism. More general population units, such as ant colonies and bird flocks, 
can also be employed and the resulting models do not differ technically from 
models treating individual organisms. 


Since individuals usually move from one location to another, it is com- 
mon in individual-based modeling to represent space explicitly. There are 
four different frameworks in which spatially-explicit individual-based models 
can be defined [5]. They differ in the way space and time are modeled: 
each can be treated either discretely or continuously. The four resulting 
frameworks have been widely studied in population ecology and they are 
considered to complement as opposed to compete with each other. In fact, 
in many cases, models of more than one types can be constructed for a sys- 
tem in question, with each approach offering a distinct perspective towards 
the understanding of the system. 


In this work, our aim is to introduce a process-calculus framework to 
enable spatially-explicit modeling of ecological systems. Indeed, process cal- 
culi also known as process algebras, first proposed in [31, 25] to aid the un- 
derstanding and reasoning about communication and concurrency, provide 
a number of features that make them suitable for capturing and reasoning 
about biological as well as ecological systems. To begin with, process calculi 
enable the construction of models by composing together smaller interact- 
ing components. This feature of modularity is especially suited towards 
the individual-based approach of modeling populations, as it enables one to 
describe the evolution of each individual of one or more populations as a 
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process and, subsequently, to compose a set of individuals (as well as their en- 
vironment) via the parallel operator into a complete ecological system. The 
semantics of the process calculus then generates automatically the complete 
state space of the system by taking into account all possible interactions and 
evolutions of the individuals/components of which the systems is composed. 
Furthermore, process calculi focus on the notion of observation, which can 
be varied according to the needs of the problem under investigation. A 
modeler can build a system incrementally and introduce less or more detail 
depending on the desired level of abstraction while, via the use of restriction 
operators, the modeler may determine the parts of the system to be stud- 
ied. Finally, the development of process calculi has been investigated into 
a wide range of directions and accompanied by a variety of analysis tech- 
niques. On the one hand, features such as time, probability and stochastic 
behavior have been extensively studied in the context of process calculi and 
they can be exploited to provide more accurate models. On the other hand, 
associated analysis techniques and tools have been developed and they can 
be used to analyze and predict system behavior. Such techniques include 
analysis of properties that can be expressed as a temporal logic formula, 
the use of equivalences for comparing the behaviors of different systems or 
for minimizing the state space of a system and, more recently, in the do- 
main of biological systems, techniques for deriving mean-field equations that 
describe the population level dynamics of individual-based models [30]. 


The process-calculus framework we propose follows the individual-based 
modeling and, in particular, it falls in the discrete-time, discrete-space class 
of Berec’s taxonomy [5]. Our process calculus, PALPS, associates processes 
with information about their location and their species. The habitat is 
defined as a graph consisting of a set of locations and a neighborhood re- 
lation. Movement of located processes is then modeled as the change in 
the location of a process with the restriction that the originating and the 
destination locations are neighboring locations. In addition to moving be- 
tween locations, located processes may communicate with each other by 
exchanging messages upon channels. Communication may take place only 
between processes which reside at the same location. Furthermore, PALPS 
may model probabilistic events, with the aid of a probabilistic choice opera- 
tor, and uses a discrete treatment of time. Finally, in PALPS, each location 
may be associated with a set of attributes capturing relevant information 
such as the capacity or the quality of the location. These attributes form the 
basis of a set of expressions that refer to the state of the environment and 
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are employed within models to enable the enunciation of location-dependent 
behavior. 

The operational semantics of our calculus is given in terms of labeled 
transition systems which combine probabilistic and nondeterministic behav- 
ior. Analysis techniques for this type of models have been widely studied in 
the formal methods literature. To allow such analysis and take advantage 
of existing tools, we describe a translation of (a subset of) PALPS into the 
PRISM language. This encoding can be employed for both simulating and 
for model checking PALPS systems using the PRISM tool [1], a model checker 
for probabilistic systems. 

In the remainder of this section, we briefly describe discrete-time discrete- 
space individual-based modeling as discussed in [5]. We continue with an 
overview of related work on the formal modeling of ecological and biological 
systems and we summarize the contributions of this article. 


Discrete-time, discrete-space, individual-based modeling. 

Individual-based modeling is an approach to model ecological systems 
which demands description of the environment and each individual living 
in it, together with its individual-individual and individual-environment in- 
teractions. It assumes that space is divided into discrete locations and that 
time runs in discrete steps. The main processes that determine the fate of 
an individual are mortality, reproduction, movement and predation. In the 
discrete-time approach, engagement in these processes is determined by a 
probability. 

At almost every step of constructing a model, at least a few alternatives 
are plausible often involving tradeoffs between making the models more 
realistic though more complex to simulate and analyze. We discuss here two 
of the parameters that involve the treatment of space and the interleaving 
of actions. 

Regarding space, discrete-space environments require a topology and 
regular lattices made up of squares are by far the most common choice. The 
size of the habitat can be considered as infinite, a choice that simplifies the 
derivation of many mathematical results, or finite, which is the choice taken 
in most simulations. In finite habitats one needs to specify the boundary 
conditions, that is, how to deal with the situation when an individual reaches 
the boundary of the habitat. There are three boundary conditions defined 
in the literature: absorbing, reflecting and periodic. Absorbing boundary 
conditions do not deal with boundary and edge effects, but focus on pop- 
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ulations in a central portion of the environment. In reflecting boundary 
conditions, individuals that hit the boundary ricochet back into the habitat 
in arandom direction. This seems realistic for animals living on an island or 
water animals living in a pond, but unrealistic for, e.g., plant seeds dispersal. 
Finally, in periodic boundary conditions, the opposite edges of the habitat 
are connected together. The argument towards the adoption of these con- 
ditions is that as one individual exits one side of the habitat another may 
enter from another one. 

Finally, another issue that is relevant to modeling ecosystems is that 
of process ordering inside each time unit. In particular, simulations carried 
out by ecologists impose an order on the events that may take place within 
a model. For instance, if we consider mortality and reproduction within 
a single-species model three cases exist: concurrent ordering, reproduction 
preceding mortality and reproduction following mortality. In concurrent or- 
dering, individuals may reproduce and die simultaneously. For reproduction 
preceding mortality, the population first reproduces, then all individuals, 
including new offspring, are exposed to death. For reproduction following 
mortality, individuals are first exposed to death and, subsequently, surviv- 
ing individuals are able to reproduce. Ordering can have significant implica- 
tions on the simulation. Thus, alternatives must be carefully studied before 
conclusions are drawn. 


Formal frameworks. Various formalisms have been proposed in the liter- 
ature for modeling biological and ecological systems. Similarly to ecosystem 
modeling, these approaches differ in their treatment of time and space and 
can be considered as supplements as opposed to rivals of each other as each 
offers a distinct view and techniques for analyzing systems. 

One strand of these formalisms is based, like PALPS, on process calculi, 
and constitutes extensions of calculi such as CCs [31], the z-calculus [32] and 
csp [25]. wsccs [45] is one of the first proposals made towards modeling 
ecosystems. It is a probabilistic extension of Ccs [31] with synchronous 
communication that has been employed in various ecological studies by the 
author and others [43, 29]. It follows the discrete-time approach to modeling 
but does not include the notion of space. WSCCS was also given a semantics 
in terms of mean field equations in [30]. This semantics captures the average 
behavior of a system over time, without computing the entire state space, 
therefore, avoiding the state-space explosion problem. 

As far as continuous time is concerned, there are various proposals in 
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the literature on stochastic process calculi. These include stochastic CCs [17], 
the stochastic pi-calculus [11] and the stochastic process algebra PEPA [24]. 
An extension of the latter, BIO-PEPA [16], has been developed for modeling 
and analyzing biochemical networks and includes features such as compart- 
ments, support for stoichiometry and general kinetic laws. 


Moving on to space, various process calculi have been proposed in the 
literature to model space in order to support the reasoning of biological 
processes. One of the most common approaches involves the notion of a 
compartment. Compartments are closed areas which may contain elements 
and other compartments and they can be used to model biological mem- 
branes or simply to represent different abstract locations. Such proposals 
include the calculus of wrapped compartments (CWC) [9], BioAmbients [41] 
as well as Brane calculi [10]. 


A different approach towards modeling biological and ecological sys- 
tems is that of P systems [40]. P systems were conceived as a class of dis- 
tributed and parallel computing inspired by the compartmental structure 
and the functioning of living cells. P-systems have been extended in various 
directions such as probabilistic P-Systems [38], continuous P-systems [37], 
dynamical probabilistic P systems [38], spatial P systems [35], and multi- 
environment P systems [13]. They have been applied to a wide range of 
applications in the field of ecology, e.g. [6, 12, 7]. 


Stochastic simulation for P systems (and other formalisms) has also 
been actively investigated. Research efforts have been based on the Gille- 
spie algorithm [22] which is currently used as the reference procedure for 
performing stochastic and discrete simulations of various biological systems. 
Of great interest is also the approximate algorithm Gillespie introduced 
in [23]. This consists of an approximate simulation strategy, referred to as 
the tau leaping method, in which, by using Poisson random numbers, it 
is possible to leap over many reaction events in a way that approximates 
the exact stochastic simulation. This strategy has also been extended in 
the spatial realm. Extensions are considered in [14] and evaluated in [26], 
whereas in [42] they were used as a basis for introducing a rewriting strategy 
for stochastic P systems of dynamically changing nested compartments and 
used for simulating ecological systems. 


Finally, we mention the calculus of looping sequences [4], and its spatial 
extension [3] and cellular automata [19, 15]. 
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Contributions. In this paper we present the process calculus PALPS, a 
process calculus for modeling spatially-explicit individual-based systems 
which adopts a discrete-time and discrete-space approach. This framework 
allows the specification of an ecosystem of multiple species residing on a 
finite topology. We model the habitat as a graph consisting of a set of lo- 
cations and a neighborhood relation that allows the modeler to define the 
intended boundary conditions. We include a special movement action that 
allows movement of individuals along the graph. 


Individuals are modeled as processes in PALPS possessing a species and 
a location that may change dynamically as computation proceeds. Individ- 
uals may engage in any of the basic processes of reproduction, dispersal, 
predation or death and they may communicate with other processes which 
reside at the same location. Probabilistic decisions are implemented through 
a probabilistic choice operator. 


Technically, PALPS can be considered as an extension of CCS with prob- 
abilistic choice, locations and location attributes. In particular, it inherits 
from ccs the notion of two-way communication according to which pro- 
cesses may interact with each other whenever two components are able to 
execute complementary actions. We believe that this type of communica- 
tion is appropriate for modeling various aspects of population behavior such 
as predation, infection and feeding: for each of these there exist two com- 
ponents able to engage in complementary behavior (e.g. predator vs. prey). 
The level of abstraction implemented by this binary (as opposed to syn- 
chronous) communication allows us to observe each local interaction and, 
for instance, to compute the number of births/infections during a time unit. 
To make this possible, interactions in PALPS, also referred to as internal 
actions in process calculi such as CCS, are labeled by the type of the action 
and the location where the action has taken place. 


As far as space is concerned, PALPS shares a similar treatment of lo- 
cations with process algebras developed for reasoning about mobile ad hoc 
networks such as [20, 27]. A significant factor that distinguishes PALPS 
from these and other process calculus approaches, is its ability to define 
location-dependent behavior. This is achieved by associating locations with 
attributes that capture information relevant to the location and form the 
basis of a set of expressions which may determine the evolution of indi- 
viduals. Examples of such attributes are the temperature, the humidity, 
and the number of preys and predators on a location. Based on such at- 
tributes, processes can determine their behavior, e.g., whether to reproduce 
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or to disperse. We illustrate the expressiveness of PALPS by constructing 
spatially-explicit individual-based models of metapopulation dispersal. 


Another contribution of our work regards the development of a method- 
ology for translating PALPS models into the PRISM language, with the prospect 
of making more advanced analysis of ecological models as opposed to just 
simulations. Most research on ecological systems models to date, focuses 
on the development of models of different metapopulation systems. One 
key question regarding these models is what one can do with the models 
other than just simulate trajectories. A possible answer to this question 
is to use model checking tools for automatically analyzing various proper- 
ties of the model. As an example, the probabilistic model checker PRISM, 
[1], has been used to investigate continuous-time P systems. Various other 
formalisms have been translated to PRISM (e.g. [34]). However, as far as 
we know, there has been no work on model checking of spatially-explicit 
individual-based models. A contribution of this paper is to make this pos- 
sible. Specifically, we propose a translation of a subset of PALPS into the 
Markov-decision-process component of the PRISM language. This subset of 
PALPS is the fragment of the language where the replication operator is 
replaced by a restricted replication operator. Essentially, this restricted op- 
erator places a bound on the number of new individuals than can be created 
of each species. This is necessary due to the fact that PRISM does not allow 
for the dynamic creation of modules, thus, the maximum total number of all 
modules/individuals must be created a-priori. We prove the correctness of 
our translation and we apply it to perform analysis of some simple examples. 


Structure of the paper. The structure of the remainder of this paper is 
as follows. In Section 2 we present the syntax and the semantics of PALPS. 
We continue to illustrate the expressiveness of the calculus by providing 
models of metapopulation dispersal in Section 3. In Section 4 we present 
a translation of PALPS into the Markov-decision-process component of the 
PRISM language. We establish the correctness of the translation and we 
overview the types of analysis that this translation makes possible on PALPS 
models. We then apply these techniques on simple examples and we explore 
the potential of the approach in Section 5. Finally, in Section 6, we conclude 
with a discussion of future work. 
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2 The Process Calculus 


In our calculus, Process Algebra with Locations for Population Systems 
(PALPS), we consider a system as a set of individuals operating in space, 
each belonging to a certain species and inhabiting a location. This location 
may be associated with attributes which describe characteristics of the loca- 
tion and can be used to define location-dependent behavior of individuals. 
Furthermore, individuals who reside at the same location may interact with 
each other by communicating upon channels, e.g. for predation, or they 
may migrate to a new location where they may continue their computa- 
tion. PALPS may model probabilistic events with the aid of a probabilistic 
operator and uses a discrete treatment of time. 


2.1 The Syntax 


In this section we formalize the syntax of PALPS which is built based on the 
following basic entities: 


S: a set of species ranged over by s, s’. 

e Loc: a set of locations ranged over by @, ¢’. The habitat of a system 
is then implemented via a relation Nb, where (¢,¢’) € Nb exactly 
when locations @ and ¢’ are neighbors. For convenience, we use Nb as 
a function and write Nb(¢) for the set of all neighbors of @. 

e Ch: a set of channels ranged over by lower-case strings. 

W: aset of attributes, ranged over by w, ~’. We write wy, for the value 

of attribute w at location @. 


Species and locations are characteristics associated with every individ- 
ual in a PALPS system. The species characteristic is static whereas the 
location characteristic is dynamic: as computation proceeds, an individual 
may change its location from ¢ to ¢’ with the restriction that (¢, 0’) € Nb. In 
turn, attributes are characteristics associated with locations and they may 
capture information such as the capacity, the temperature or the quality of 
the location. They form the basis of the set of expressions of the language 
which is defined below. 


Expressions. PALPS employs two sets of expressions: logical expressions, 
ranged over by e, and arithmetic expressions, ranged over by w. These 
expressions are intended to capture environmental situations which may 
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affect the behavior of individuals. Expressions e and w are constructed as 
follows: 


e s= true | n7e |e, Ae | wc 
wos= c| vQ@& | sQe* | @ | op,(w) | ope(wi, we) 


where c is a real number, me {=,<, >} and & € Loc U {myloc}. 

To begin with, logical expressions e are built using the propositional 
calculus connectives, as well as comparisons between an arithmetic expres- 
sion w and a constant c (e.g., s;}@é@ + s2Q@é@ > 1). Arithmetic expressions 
include three special expressions interpreted as follows: Expression ~@é* is 
equal to the value of attribute ~ at location @*. Expression (s@é*) is equal 
to the number of individuals of species s at location @*, and expression @é* 
denotes the total number of individuals of all species at location £*. 

Location €* can be an arbitrary location or the special location myloc. 
This latter label is employed to bestow individuals with the ability to express 
conditions on the status of their current location no matter where that might 
be, as computation proceeds. Specifically, myloc refers to the actual location 
of the individual in which the expression appears and it is instantiated to 
this location when the condition needs to be evaluated (see rule (Cond) in 
Table 3). In conclusion, arithmetic expressions are the set of all expressions 
formed by arbitrary constants c, quantities ~@0*, s@é*, @é*, and the usual 
unary and binary arithmetic operations (op; and op.) on real numbers. 
Logical expressions and arithmetic expressions are evaluated within a system 
environment (as defined in Tables 1 and 2). 


Processes. The syntax of PALPS is given at three levels: (1) the individual 
level ranged over by P, (2) the species level ranged over by R, and (3) the 
system level ranged over by S. Their syntax is defined via the following 
BNFs: 


Po 0) eres | >opiPi | cond (6) > Piys.c5@n & P,) | -¢ 
tel wel 

Ro o:= !rep.P 

S := 0 | Pxs,é) | Rvs) | Si]So | S\L 


where L C Ch, J is an index set, p; € (0,1) with S0,-;p; = 1, €1,..-,€n, 
is a set of logical expressions such that e; V...V en, = true, C’ ranges over 
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a set of process constants C, each with an associated definition of the form 


C= P, and 


nese |-a)|-gat |f. 


Beginning with the individual level, P can be one of the following: 
Process O represents the inactive individual, that is, an individual who has 
ceased to exist. Process Se 1m-P; describes the nondeterministic choice be- 
tween a set of action-prefixed processes: it can execute any of the activities 
7, and proceed as the respective P;. We write 7,.P, + 72.P2 to denote the 
binary form of this operator. In turn, an activity 7 can be an input action 
on a channel a, written simply as a, a complementary output action on a 
channel a, written as G, a movement action with destination ¢, go@, or the 
time-passing action, written as \/. Actions of the form a, and @, a € Ch, are 
used to model arbitrary activities performed by an individual; for instance, 
eating, predation and reproduction. The tick action \/ measures a tick on a 
global clock and is used to separate the rounds of an individual’s behavior. 
Essentially, the intention is that in any given time unit all individuals per- 
form their available actions possibly synchronizing as necessary until they 
synchronize on their next \/ action and proceed to their next round. 

Process »),< ;pi:P; represents the probabilistic choice between processes 
P;, i © I. The process randomly selects an index 7 € J with probability p;, 
and then evolves to process P;. We write p,:P, @p2:P», for the binary form of 
this operator. The conditional process cond (e1 > P1,...,@n & Pn) presents 
the conditional choice between a set of processes: it behaves as P;, where i 
is the smallest integer for which e; evaluates to true. Note that this choice is 
deterministic. Finally, process constants provide a mechanism for including 
recursion in the calculus. 

Moving on to the species level, we employ the special species process 
R defined as !rep.P. This process is a replicated process which may always 
receive input through channel rep and create new instances of process P, 
where P is a new individual of species R. Such inputs will be provided by 
individuals in the phase of reproduction via the complementary action rep. 

Finally, population systems are built by composing in parallel located 
individuals and species. An individual is defined as P:(s, €), where s and ¢ 
are the species and the location of the individual, respectively. A species 
is given by R:(s), where s is the name of the species. Finally, S\Z models 
the restriction of the use of channels in set LZ within S. As a syntactic 
shorthand, we write P:(s,@,n) for the parallel composition of n copies of 
process P:(s, £). 


130 A. Philippou, M. Toro, M. Antonaki 


Example 1. We consider a simplification of the model presented in [44] 
which studies the reproduction of the parasitic Varroa mite. This mite usu- 
ally attacks honey bees and it has a pronounced impact on the beekeeping 
industry. In this system, a set of individuals reside on an n x n lattice of 
resource sites and go through phases of reproduction and dispersal. Specif- 
ically, the studied model considers a population where individuals disperse 
in space while competing for a location site during their reproduction phase. 
They produce offspring only if they have exclusive use of a location. After 
reproduction the offspring disperse and continue indefinitely with the same 
behavior. In PALPS, we may model the described species s as R 1 -ep.P, 
where 


1 
KR & > 1: gol.,/.cond (sQmyloc = 1 > Pi; true \/.Po) 
LENb(myloc) 


PR & p:rep.,/.Po © (1 — p):rep.rep./. Po 


According to the definition, during the dispersal phase, an individual 
moves to a neighboring location which is chosen probabilistically among the 
four neighboring locations on the lattice of the individual. Subsequently, 
the conditional construct allows to determine the exclusive use of a location 
by an individual. The special label myloc is used to denote the actual 
location of an individual once the individual is instantiated within a system 
definition. In case of exclusive usage, the flow of control process according 
to P, which models the probabilistic production of one or two children of 
the species. A system that contains two individuals at a location 2 and one 
at location ¢’ can be modeled as 


System © (Po:(s, £, 2)|Po:(s, £’)|R:(s))\{rep}. 


Let us now extend the example into a two-species system. In particular, 
consider a competing species s’ of the Varroa mite, such as the pseudo- 
scorpion, which preys on s. To model this, we may define the process 


def 
R' ='lrep!.Qo, where 


Qo = prey./.Qit V.Q2 
Qi = rep V/.Qo 
Qo S prey./.Qit+ V0 


Simulation and Verification in a Process Calculus for 
Spatially-Explicit Ecological Models 131 


An individual of species s’ looks for a prey. If it succeeds in locating one, 
then it produces an offspring. If it fails for two consecutive time units it 
dies. 

To implement the possibility of preying on the side of species s, the 
definition must be extended by introducing the complementary input actions 
on channel prey at the appropriate places: 


il 
Po e ~ 1: (go l../.cond (s@myloc = 1 > P;; true > ./.Po) + prey.0) 
LENb(myloc) 
P, & p:(rep../.Po + prey.0) & (1 — p):(rep.rep.,/.Po + prey.0) 


2.2 The Semantics 


The semantics of PALPS is defined in terms of a structural operational se- 
mantics given at the level of configurations of the form (£, 5), where E is an 
environment and S is a population system. The environment F is an entity 
of the form E Cc Loc x 8 x N, where each pair @ and s is represented in E at 
most once and where (¢, s,m) € FE denotes the existence of m individuals of 
species s at location ¢. The environment F plays a central role in evaluating 
expressions. 

The satisfaction relation for logical expressions — is defined inductively 
on the structure of a logical expression, as shown in Table 1. It depends 
on the evaluation function for arithmetic expressions val (£,w) defined in 
Table 2. 


Table 1: The satisfaction relations for logical and arithmetic ex- 
pressions 


E/E true always 
E |= -7e if and only if -=(E Fe) 
EFeAeg ifandonlyif BEEeAE Ee: 


EEwme  ifandonlyif val(£,w)™e 


Before we proceed to the semantics we define some additional opera- 
tions on environments that we will use in the sequel: 


Definition 1. Consider an environment EF, a location @ and a species s. 
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Table 2: The evaluation relation for arithmetic expressions 


val (£7, c) = <¢ 

val (E, @e) = we 

val (EL, sQ£) = n,(f,s,n)EeE 

val (EL, sQ£) = 0,(¢s,n)¢E 

val (EL, @£) = Dees ts, (4,8, %s) € E 

val (7, op; (w)) =  op,(val (E,w)) 

val (Z,ops(wi,we)) = op,(val (EL, wz), val (E,we)) 


e E'&(s,£) increases the count of individuals of species s at location @ 
in environment E by 1: 


_f BU{(és,m+1)} if B= E'U{(é,s,m)} for some m 
Beh { EU {(é,s,1)} otherwise 


e £6 (s, 2) decreases the count of individuals of species s at location ¢ 
in environment FE by 1: 


E’U{(é,s,m—-1)} if B= E’U{(e,s,m)},m>1 
Eo(s,4)=<¢ E' if FE = FE’ U{(é,s,1)} 
L otherwise 


We may now define the semantics of PALPS, presented in Tables 3 
and 4, and given in terms of two transition relations: the non-deterministic 
relation —+, and the probabilistic relation —>,. A transition of the form 
(E,S) “+, (B', 8") means that a configuration (E,S) may execute action 
and become (E£’, $’). A transition of the form (E,S) +, (E’, 5’) means 
that a configuration (EZ, S) may evolve into configuration (E’, S’) with prob- 
ability w. Whenever the type of the transition is irrelevant to the con- 
text, we write (E,S) > (E’,S’) to denote either (E,S) +, (E’,S’) or 
(E,S) +, (E’,S’). Action js appearing in the non-deterministic relation 
may have one of the following forms: 


e aQ@é and @@é denote the execution of a and @ respectively at location 
é; 

e Taae denotes an internal action that has taken place on channel a at 
location £. This may arise when two complementary actions take place 
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at the same location @ or when a move action take place. Note that 
recording the location and the channel of a synchronization in internal 
actions, is a point on which PALPS departs from other process calculi 
featuring binary communication. The intention behind this design 
decision is to enable to observe and reason about behavior of a system 
pertaining to synchronizations. 

e ,/ denotes the time passing action. 


Rules for individuals. The rules of Table 3 prescribe the semantics of 
located individuals in isolation. The first four rules define non-deterministic 
transitions. The fifth axiom defines a probabilistic transition, and the last 
two rules refer to both the non-deterministic and the probabilistic case. All 
rules are concerned with the evolution of the individual in question and 
the effect of this evolution to the system’s environment. A key issue in 
the definition of the rules is to preserve the compatibility of P and FE as 
transitions are executed. We consider each rule separately: 


e Axiom (Nil) specifies that the 0 process may execute the time consum- 
ing action ,/. This axiom allows for time-progress in a system with 
inactive individuals. 

e Axiom (Tick) specifies that a \/-prefixed process will execute the time 
consuming action \/ and then proceed as P. The state of the new 
environment depends on the state of P. If P = 0 then the individual 
has terminated its computation and it is removed from E (see the 
definition of E”*-"), If P 4 0 then E remains unchanged. 

e Axiom (Act) specifies that 7.P executes action 7@¢ and evolves to P. 
Note that the action is decorated by the location of the individual 
executing the transition to enable synchronization of the action with 
complementary actions taking place at the same location (see rule 
(Par2), Table 4). This axiom excludes the cases of 7 = gof and n = \/ 
which are treated separately. 

e According to Axiom (Go), an individual may change its location. This 
gives rise to action Tg,ag and has the expected effect on the environ- 
ment F. 

e Rule (NSum) describes the behavior of a nondeterministic choice: any 
of the available summands may be selected and executed. 

e Rule (PSum) expresses the semantics of probabilistic choice: a process 
is chosen probabilistically leading to the appropriate continuation. If 
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Table 3: Transition rules for individuals 


Nil)  (B,0:(s, €)) “sy, (E, 0:(s, £)) 

Tick) (BE, /.P:(s, 0)) “on (BP, P:(s, £)) 

Act) (BE, n.P:(s, €)) 7S, (BPS, P:(s, £)) n# gol, J 
2h (EP SEL Pls 0) (¢,/) €Nb 


= 
£)) +n (E", Pir(s,£)) 
(E, Si.Pi:(s, 0) Sn (B', Pi:(s, 2) 


PSum) (EB, WierPitPi(s, £)) Srp (BP **, Pi:(8,) 
(EB, P:(s, £)) > (E", P's, £)) cas 
(E,C:{s, £)) <> (E', P's, £')) 

(E, Pi:(s, €)) —> (E', Pi:(s,€)), ERet6E eth j <i 
(E. cond (€; > Pi,...,€n & Pn):(s,£)) > (E", Pi:(s, @)) 


P:(s, €) 


Ee(s,0) ifP=0 
: Pysjk _ ) 
where ~ { E otherwise 


EP tl’ — (Ee (s,0)) @(s, 0)) Ps" 


the resulting state of the individual, namely P;, is equal to 0, then the 
individual is removed from the environment EF. 

e Rule (Const) expresses the semantics of process constants. 

e Finally, rule (Cond) stipulates that a conditional process may perform 
an action of continuation P; assuming that e;| @ evaluates to true and 
all e; | 2, 7 < 7 evaluate to false. Note that we write e | ¢ for the 
expression e with all occurrences of myloc substituted by location £. 


Rules for systems. We may now move on to Table 4 which defines the 
semantics of system-level operators. The first two rules define the semantics 
for the replication operator, the next five rules define the semantics of the 
parallel composition operator, and the last rule deals with the restriction 
operator. 

According to axioms (R_Tick) and (R_Rep), a species process may idle or 
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Table 4: Transition rules for systems 


(R-Tick) R =!rep.P:(s) 
(E, R) +n (ER) 
(R_Rep) R =!rep.P:(s), € € Loc 


(E, R) "2S", (£ @ (s, 8), P:ls, ®|R) 


(Party (E251) “en (2', 81), (E, S2) Ap AV 
(EB, S1|S2) sa (E’, S}|S2) 


(Par2) (2 51) on (E181), (E,S2) nu (Ex, $4) 
(E, S1|S2) =%,, (E ® (£1, E2), $1 |S) 
(Par3) (E, $1) ~» (E1, 51), (E, S2) *+y (E2, 95) 
(E, S1|S2) “**?, (E @ (Ei, E2), 51|$3) 
(Par4) (E, 51) rp (E£", 51), (E, S2)p 
(E, S\|S2) +p (E’, $1 |S2) 
(rime. “8a “sn (Er, $4), (EB, S2) “on (Ea, 95) 
(E, $1|S2) “+n (E ® (Ei, Ea), 54]94) 
(Res) (E, S) > (E',S"),a ¢ {a@e, a@éla € L} 


(B,S\L) SS (EUS )\\L 


it may engage in action repQ@é for any location ¢ and create a new individual 
P of species s at location @. 

Rules (Parl) - (Par4) specify the semantics of parallel composition. 
Note that the symmetric versions of these rules are omitted. According 
to (Parl), if a component may execute a non-deterministic transition and 
no probabilistic transition is enabled by the other component (denoted by 
(E, S2)p), then the transition may take place. If the parallel components 
may execute complementary actions on some channel a and some location 
é, then they may synchronize with each other producing action T,@¢ (rule 
(Par2)). If both components may execute probabilistic transitions then they 
may proceed together with probability the product of the two distinct proba- 
bilities (rule (Par3)). If exactly one of them enables a probabilistic transition 
then this transition takes precedence over any non-deterministic transitions 
of the other component (rule (Par4)). 
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Note that in case that the components proceed simultaneously then 
the environment of the resulting configuration should take into account the 
changes applied in both of the constituent transitions (rules (Par2), (Par3) 
and (Time)). This is implemented by FE ® (EF), E2) as follows: 


E® (Fy, E2) = {(2,s,m+%4 + ig) | (¢,s,m) € EB, 
(¢,s,m + 74) € Fy, (¢,s,m + i2) € E9,%1,12 € Z} 


Rule (Time) defines that parallel processes must synchronize on ,/ ac- 
tions. This allows one tick of time to pass and all processes to proceed to 
their next round. Finally, rule (Res) defines the semantics of the restriction 
operator in the usual way. 


Initial configuration. Based on this machinery, the semantics of a sys- 
tem S is obtained by applying the semantic rules to the initial configuration. 
The initial configuration, (£,S'), is such that (€,s,m) € E if and only if S 
contains exactly m individuals of species s located at ¢ of the form P:(s, £), 
where P # 0. In general, we say that E is compatible with S whenever 
(¢,s,m) € E if and only if S contains exactly m individuals of species s 
located at @. It is possible to prove that the defined semantics preserves 
compatibility of configurations [2]: 


Lemma 1. Whenever (E,9) + (E',S’) and E is compatible with S, then 
E' is also compatible with S’. 


3 Examples 


During the last few decades, the theory of metapopulations has been an 
active field of research in Ecology and it has been extensively studied by 
conservation scientists and landscape ecologists to analyze the behavior of 
interacting populations and to determine how system parameters may in- 
fluence various aspects of these systems such as local and global popula- 
tion persistence and species evolution. Population dispersal is one such 
phenomenon that has been of special interest to ecologists. It affects the 
long-term persistence of populations, the coexistence of species and genetic 
differentiation between subpopulations and understanding this process is 
essential for obtaining a good understanding of the behavior of metapopula- 
tions. The evolution of dispersal has received much attention by scientists 
and it has been studied in connection to various parameters such as the 
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Figure 1: The sequence of events in the lifetime of a dispersing species 


connectivity of the habitat on which a metapopulation exists, patch quality 
and local dynamics. 

In this section, we describe two examples relating to metapopulation 
dispersal through which we illustrate how our calculus can be used to con- 
struct models of this phenomenon. 


Example 2. This example is motivated by the spatially-explicit individual- 
based model of [46]. In this work the authors construct a fairly simple 
model of metapopulation dispersal which departs from previous works in 
that, unlike previous models of metapopulation dispersal which tended to be 
deterministic and at the level of population densities, the model constructed 
is stochastic and individual-based. 

According to this study, a set of genotypes co-exist within a habitat and 
differ only in their propensity to disperse. The metapopulation is composed 
of n x n subpopulations inhabiting a set of patches arranged on a square 
lattice with cyclic boundaries, so that individuals leaving the “top” or “right- 
side” of the world reappear on the “bottom” or “left-side” respectively and 
vice versa. Each patch is associated with a so-called patch quality related 
to the capacity of the patch. 
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The behavior of an individual of the genotypes under study is illus- 
trated diagrammatically in Figure 1. To begin with, an adult individual 
produces A descendants. Subsequently, a phase of competition takes place 
between the juveniles of the population of which a fraction survives. Each 
surviving offspring may disperse according to a probability of dispersal dis- 
tinct to its genotype. In case it disperses, the neighboring patch it moves 
to is selected with equal probability among all neighbors. We point out 
that the percentage of offspring surviving juvenile competition at patch @ is 
given by ye = (1 + ag- Ne)’, where ay is the measure of the patch quality, 
N¢ is the number of individuals residing at patch @ and ( is a constant that 
relates to the degree of competition. 

This metapopulation can be modeled in PALPS as follows. We consider 
the set of locations (i, 7), 1 < 7,7 <n, where two locations (7,7) and (k,l) 
are neighbors if they are adjacent on the grid. Further, we use the location 
attribute ag as a measure of the quality of the patch at @. Then, genotype 7 
with some constant probability of dispersal p; and A = 3 can be defined as 
the species process R; =!rep;.J;, where 


A de ep, .Tep,.Fep;.0 Adult Individual 

x = qiiS; ® (1 — q):0 Juvenile 

S; def pj:D; © (1 — pj) fA; Surviving Juvenile 
Di def See : go b.\/ Aj Dispersing Juvenile 


and q; the probability of survival of juvenile competition is given by q; = 
(1 + ag- @0)°. Then a system can be modeled as the composition of the 
various genotypes as well as the individuals of the initial population under 
study: 


System © ((Ry:(1) | ..-|Retk)| [Ais (é,,ma))\{repi,..repe}. 


1<i<k,leELoc 


Analysis in this model may focus on the effect that the dispersal rates, 
the degree of competition and/or the patch quality may have on population 
dispersal. 


Example 3. As another more complex example, let us consider a model of 
wood thrush dispersal, initially proposed in [47] and expanded upon in [33]. 
This model considers three types of birds: adult breeders, adult floaters, 
and juveniles which are birds in their first year of life. 
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According to this model, adult breeders produce an offspring at a rate 
dictated by various system parameters such as clutch size, nest predation 
and parasitism rates which we denote as ry. Following reproduction, each 
individual has a probability of dying before the next time step which is 
higher in juveniles and adult floaters in comparison to adult breeders. We 
write q, q; and q, for the mortality rates of breeders, juveniles and floaters, 
respectively. If following mortality a habitat patch has more birds than 
its capacity allows, then dispersal will occur according to a probability de- 
termined by the size of the patch and the distance between neighboring 
patches. This probability is higher in floaters and juveniles in comparison 
to adult breeders who exhibit a high site fidelity. We write p,, p; and pf 
for the dispersion rates of breeders, juveniles and floaters, respectively. If 
a bird reaches a patch with available capacity then it will settle. If not, 
then it will either attempt to disperse to another patch or it will become a 
floater depending on whether it has reached its maximum number of disper- 
sal events. Once dispersal has occurred, the juveniles become adults and 
the model begins another cycle. This sequence of events in the behavior of 
the populations is presented diagrammatically in Figure 2. 


This metapopulation can be modeled in PALPS as follows. We consider 
the set of locations Loc and an associated predefined neighbor function. 
We also assume the existence of a set of probabilities {pj,;}i;eLoc where pi,j 
represents the probability of dispersal from patch 7 to patch 7. Finally, we 
introduce the location attribute cp as a measure of the capacity of patch @. 
Then, the wood thrush species can be modeled by the process R =!rep. Juv, 
where the behavior of a juvenile individual J; is described by the following 
equations: 


Tug os qj: JCo ® (1 —q;) : 0 Juvenile survival 
JICo det cond (Q@myloc > Cmyloc & JDo, true > /.AB) Check patch capacity 
JDo det pj: JA1®(1—p;): /.AB Decide whether to disperse 
JA, %& Serre ss Bind" gol.JICy Dispersal attempt 1 
IC, def cond (@myloc > Cmyloc > JD1, true > /.AB) Check patch capacity 
JID, def pj: JA2 ®(1—p;): /.AB Decide whether to disperse 
cee def hae 2c) Pmyloc,e ? JO L.IC2 Dispersal attempt 2 
ef 


JC. = cond (@myloc > Cmyloc > \/-F'l, true > \/.AB) Become floater or adult 
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AB |= 1r:Tep.BS®(1—rb): BS Breeder reproduction 
BS |= q@:BCo@(1—q):0 Breeder survival 
BCo = cond (@myloc > Cmyloc > BDo, true > \/.AB) Check patch capacity 
BDo = py: BA: ®(1—py): /.AB Decide whether to disperse 
BA, & ce b(myloc)Pmyloc,£ : gol.BC, Dispersal attempt 1 
BC, = cond (Q@myloc > Cmyloc > BD1, true > \/.AB) Check patch capacity 
BD, = py: BAs@®(1—py): /.AB Decide whether to disperse 
BA, & cen (myloc)Pmyloc,£ : gol.BC2 Dispersal attempt 2 
BCz = cond (@myloc > Cmyloc > /.L'l, true > \/.AB) Floater or adult 
Fl © af: FCo®(1—qr): 0 Floater survival 
FCo = cond (@myloc > Cmyloc > F'Do, true > \/.Fl) Check patch capacity 
FDo © ps: FA: @(1—pr): V-FI Decide whether to disperse 
ay Cae ce b(myloc)Pmyloc,£ : gol. FC, Dispersal attempt 1 
FC, © cond (@myloc > cmyloc & FD1, true > /.Fl) Check patch capacity 
FD, = py: FAr@(1— py): -FI Decide whether to disperse 
FA, & yo cen (myloc)Pmyloc,e : g0l./.Fl Dispersal attempt 2 


As before, the system can be modeled as the composition of the species as 
well as the various individuals that form the study: 


System a (R:(1) | II AB:(1,£, nb) | II Juv:(1, £,n5) 
leLoc £e€Loc 


| J] Fld, 4n4))\{rep} 


LeELoc 


Varying the model parameters (such as the habitat topology and the patch 
quality) may allow an analysis of the effects of the parameters on patch and 
metapopulation persistence. 


4 ‘Translating PALPS into PRISM 


In this section we turn to the problem of model checking PALPS models. To 
begin with, we observe that the operational semantics of PALPS gives rise to 
transition systems that can be easily translated to Markov decision processes 
(MDPs). We recall that Markov decision processes are a type of transition 
systems that combine probabilistic and non-deterministic behavior. As such, 
model checking approaches that have been developed for MDPs can also be 
applied to PALPS models. 
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Given this observation, in this section we explore the possibility of 
employing the probabilistic model checker PRISM to perform analysis of the 
semantic models derived from PALPS processes. Specifically, we propose a 
translation scheme for encoding PALPS systems into PRISM models on which 
model checking can be performed. In the next section, we will illustrate 
our approach by performing analysis on the PALPS system presented in 
Example 1. 

PRISM is a probabilistic model checker for Markov decision processes, 
discrete time Markov chains, and continuous time Markoc chains. Using 
PRISM it is also possible to generate random sample paths of execution for 
simulation. For our study we are interested in the MDP support of the tool 
which supports model checking though probabilistic computational time 
logic [1]. Specifically, we will propose a method for translating PALPS sys- 
tems into the MDP subset of the PRISM language yielding models on which 
simulation and analysis can be performed via the PRISM model checker. 

There are two alternatives for generating an input for PRISM given a 
PALPS model. First, it is possible to explicitly construct the Markov decision 
process corresponding to the model and directly specifying the transition ma- 
trix to the model checker. Second, it is possible to reproduce/translate the 
model in question into the PRISM language. According to PRISM’s manual, 
the latter method is typically more efficient than the first alternative. This 
is because PRISM is a symbolic model checker and the underlying data struc- 
tures used to represent the model, function better when there is a high-level 
structure and regularity to exploit. Therefore, we chose to encode PALPS 
models into the PRISM language. In the remainder of this section, we will 
give a brief presentation of the PRISM language, present an encoding of (a 
subset of) PALPS into PRISM and prove its correctness. 


4.1 The PRISM Language 


The PRISM language is a simple, state-based language, based on guarded 
commands. A PRISM model consists of a set of modules which can interact 
with each other on shared actions following the CSP-style of communication. 
Each module possesses a set of local variables which can be written by the 
module and read by all modules. In addition, there are global variables 
which can be read and written by all modules. The behavior of a module is 
described by a set of guarded commands. When modeling Markov decision 
processes, these commands take the form: 


Simulation and Verification in a Process Calculus for 
Spatially-Explicit Ecological Models 143 


[act] guard py : uw +... + Pm iUm;3 


where act is an optional action label, guard is a predicate over the set of 
variables, p; € (0, 1] and u; are updates of the form: 


fos = Ui) &... & (x, = Uik) 


where u;,; is a function over the variables. Intuitively, such an action is 
enabled in global state s if s satisfies guard. If a command is enabled 
then it may be executed in which case, with probability p;, the update u; 
is performed by setting the value of each variable x; to ui,j(s) (where x; 
denotes the new value of variable «;). 

A model is constructed as the parallel composition of a set of modules. 
The semantics of a complete PRISM model is the parallel composition of all 
modules using the standard CSP parallel composition. This means that all 
the modules synchronize over all their common actions (i.e., labels). For 
a transition arising from synchronization between multiple processes, the 
associated probability is obtained by multiplying those of each component 
transition. Whenever, there is a choice of more than one commands, this 
choice is resolved non-deterministically. We refer the reader to [1] for the 
full description and the semantics of the PRISM language. 


4.2 Encoding PALPS into the PRISM Language 


As observed in [34], the main challenge of translating a ccs-like language 
(like PALPS) into PRISM is how to map binary Ccs-style communication over 
channels to PRISM’s multi-way (CSP style) communication. Our approach for 
dealing with this challenge, similarly to [34], is to introduce a distinct action 
for each possible binary, channel-based, communication which captures the 
channel as well as the sender/receiver pair. The two other actions in PALPS, 
namely the tick action and the movement action, can be easily handled as 
follows: In the case of the tick action, the encoding is straightforward as 
its semantical treatment in PALPS is via synchronous communication which 
can be represented directly in PRISM. In the case of movement actions, 
it is possible to introduce a PRISM command that performs the necessary 
updates on the state of the model. 

To translate PALPS into the PRISM language, we translate each process 
into a module. The execution flow of a process is captured with the use of a 
local variable within the module whose value is updated in every command 
in such a way as to guide the computation through the states of the process. 
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Then, each possible construct of PALPS is modeled via a set of commands. 
For example, the probabilistic summation is represented by encoding the 
probabilistic choices into a PRISM guarded command. Non-deterministic 
choices are encoded by a set of simultaneously enabled guarded commands 
that capture all nondeterministic alternatives, whereas the conditional state- 
ment is modeled as a set of guarded commands where the guard of each 
command is determined by the expressions of the conditional process. 

Unfortunately, the replication operator cannot be directly encoded into 
PRISM since the PRISM language does not support the dynamic creation of 
modules. To overcome this problem, we consider a bounded replication 
construct of the form !”P in which we specify the maximum number of 
P’s, namely m, that the can be created during computation. We note that, 
in practice, the value of m can be selected by estimating a bound of the 
maximum size of the population, or it can be determined by the size of the 
state-space of the resulting model. 

We now consider the main ideas of translating PALPS into the PRISM 
language via an example. 


Example 4. Consider a habitat consisting of four patches {1, 2,3, 4}, where 
Nb is the symmetric closure of the set {(1,2), (1,3), (2,4), (3,4)}. Further- 
more, suppose that each location @ has an attribute cy determining the 
capacity of the location in terms of the number of individuals it may sup- 
port. Let s be a species residing on this habitat defined according to the 
bounded replication R: 


R = !"rep.P 
P © cond (sQmyloc > c@myloc > P;, s@myloc < c@myloc > P2) 
f 1 
Pp, = ~ 5 : gol../P 
LENb(myloc) 
d 


pe Tep../.P 


According to the definition, an individual of species s begins by check- 
ing whether the location has exceeded its capacity. If so then it migrates 
with equal probability to one of the neighboring locations, otherwise, it 
produces one offspring and it returns to its initial state. Now, consider a 
system initially consisting of two individuals, one at location 1 and one 
at location 4: 


System © (P:(s, 1) | P:(s,4) | R:(s))\{rep} 
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In order to translate System in the PRISM language we first need to 
encode global information relating to the system. This consists of four 
global variables that record the initial populations of each of the locations 
and the values of attributes cy. For simplicity, let us assume that c = 2 for 
all locations. We also include a global variable 7 that measures the inactive 
individuals still available to be triggered. Initially 7 = m. 


const int c = 2; 

global s1: [0,m+2] init 1; 
global s2: [0,m+2] init 0; 
global s3: [0,m+2] init 0; 
global s4: [0,m+2] init 1; 
global i: [0,m+2] init m; 


We continue to model the two individuals P:(s,1) and P:(s,4). Each 
individuals will be described by a module. In Figure 3, we may see the 
translation of individual P:(s,1). We observe, that its species variable sp1 
is set to 1, a constant that identifies the species, its location variable locl 
is set to 1 and variable stl, recording its state, is set to 1, the initial state 
of the module. Overall, the module has 5 different states. From state 1 
two commands are enabled: one leads to state 2 and the other to state 4. 
From state 2, reproduction is implemented in two steps: in the first step 
the number of individuals of the current location is increased by one and 
stl becomes 3, and in the second step, the offspring is created by synchro- 
nization of the module and an inactive individual module via action rep_1_1. 
Note that these two moves cannot be merged into one due to the restriction 
of PRISM that commands which synchronize with other modules (such as 
rep_1i) cannot modify global variables (such as sz). On the other hand, 
from state 4, a probabilistic transition is enabled during which the module 
changes its location (variable locl) to one of the neighboring locations and 
updates the populations of the relevant locations accordingly. In both cases, 
the flow of control is determined by the assignment st = 5. Whenever st = 5, 
the module is willing to synchronize via action tick, which we assume to be 
an action shared by all modules and via which a tick on the global clock is 
measured. Then the module resumes the initial state (st = 1). 

Individual P:(s,4) may be defined similarly. Note that the module 
encoding this process is identical to module P1 with the exception that we 
rename the variables and the initial value of the location variable. 

This leaves us with the encoding of R, the component that implements 
replication of individuals. As we have already discussed, we achieve this via 
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module Pi 


sti: [1..5] init 1; 
loci: [1..4] init 1; 
const int spl = 1; 


11 (st1=1)&(si<c) -> (st1’=2); 


[] (st1=2)&(loc1=1)&(i>0) -> (s1’=s1+1)&(i’?=i-1) &(st1’=3); 
[rep_1_3] (st1i=3) -> (st1’=5); // Activate module 3 
... // All activation possibilities are enumerated 


[] (st1=2)&(1loc1=2)&(i>0) -> ... 
[] (st1=2)&(1oc1=3)&(i>0) -> ... 
[] (st1=2)&(1loc1=4)&(i>0) -> ... 


1] (st1=1)&(si>=c) -> (st1’=4); 

(] (st1=4)&(loct=1) -> 

0.5: (st1’=5)&(loc’=2)&(s1’=s1-1)&(s2’=s2+1) 

+ 0.5: (st1’=5)&(loc’=3) &(s1’=s1-1)&(s3’= s3+1); 
[] (sti=4)&(loci=2) -> ... 

[] (sti=4)&(loci=3) -> ... 

[] (sti=4)&(loci=4) -> ... 


[tick] (sti=5) -> (st1’=1); 
endmodule 


Figure 3: PRISM code for an active individual 


bounded replication which makes an assumption on the maximum number of 
new individuals that can be created in a system. Given this assumption, our 
model must be extended by an appropriate number of inactive individuals 
awaiting for a trigger via a rep_i_j action as illustrated in Figure 4. 

Thus, the inactive individual modeled by module P3, awaits to synchro- 
nize with any of the remaining modules 1,...,.maz, in which case it inherits 
the location of the synchronizing module and it sets st3 = 1 so that it may 
begin to execute the code of an active individual, presented in Figure 3, 
with the variables appropriately renamed. 


4.3. Formal Translation 


In this section, we will formalize the intuitions of the previous example into 
a formal translation of PALPS into PRISM and we will prove its correctness. 
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module P3 
st3 : [0..5] init 0; 
loc3: [1..4] init 1; 


const int sp3 = 1; 


[tick] (st3=0) -> (st3’=0); 
[rep_1_3] (st3=0) -> (st3’=1)&(1oc3’=1loc1); 


[rep_max_3] (st3=0) -> (st3’=1)&(1loc3’=loc_max) ; 


// Here we append the code of an active P individual 
// with the variables appropriately renamed 


endmodule 


Figure 4: PRISM code for an inactive individual 


Consider a PALPS model. This consists of a set of locations, L = 
{1,...,k}, a set of attributes, O = {0),...,8mn}, and a value of each at- 
tribute at each location, the neighborhood relation Nb and a process System. 
We assume, for the reasons already discussed in the previous example, that 
all replication processes are bounded and have the form R =!"rep.P, thus 
allowing the creation of up to n individuals of the specific species. We also, 
assume that all rep channels are restricted within System. Then, the PRISM 


model is constructed as follows: 


e For each species s, the model contains the & global integer variables 
$1,-.-, Sk, capturing the number of individuals of species s for each of 
the locations. The variables are appropriately initialized based on the 
definition of System. 


For each attribute 6 and each location £, the model contains a constant 
that records the values of 0,. 


Each (active) process P:(s, 0) of System becomes a PRISM module with 
a constant spp = s, a variable locp = @ and a variable stp, with range 
1,...,|P|, which records the current state of the individual and where 
|P| is the number of states that P may evolve to. The body of the 
module is the translation of process P into guarded commands. 


e Each species definition R:(s) =!"rep.P of System becomes a sequence 
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of n PRISM modules, Pyi1,..-,;Pe4n, where x is the number of indi- 
viduals of species s in the initial state of System. In our model, we 
introduce a variable i, that records the current number of inactive 
individuals of species s: 


is : [0..n] init n; 


Each inactive module Py possesses a constant spy = s, a variable 
locy which is not initialized, and another variable st,, with range 
0,...,|P|, which corresponds to the current state of the individual 
and where |P| is the number of states that P may evolve to. Note 
that st, = 0 corresponds to the state where the inactive individual is 
awaiting activation by one of the active individuals of species s. To 
capture this, we include the following commands: 


module Py 

sty : [0..|P|] init 0; 

locy : [1..m]; 

const int sy = 8; 

[tick](sty = 0) —> (st, = 0) 

[repiy](sty = 0) —> (sti, = 1)&(locy = loci); 


[repriny](sty = 0) —> (st, = 1)&(locy = loc); 


//Here we append the translation of P 


endmodule 


Thus, Py, may be activated by any of the modules P,,..., Pr4n and 
then proceed according to the translation of process P into guarded 
commands. 


We now continue to describe how a process at the individual level of 
PALPS can be translated to a sequence of PRISM commands. We denote the 
translation of a process P as | P] and we define it inductively on the structure 
of P. Note that, for convenience, we write sg to for the state (integer value) 
associated with process @. Finally, in the translations below, we assume 
that we are working within a module with identifier x and variables st, and 
loc. 
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Case 1: Q = gol.P. We translate the process by including the command 
[| (st = sg)&(loc = a) —> (st’ = sp)&(loc’ = I)& (st, = sta — 1)& (st) = st) + 1); 

for each (a,l) € Nb and then appending the translation of P. 

Case 2: Q =a.P. We translate the process by including the command 


[ay,2] (st = sq) &(loc = locy) —> (st' = sp); 


for each location / and each module y that may perform an output on 
channel a. We then append the translation of P. 
Case 3: Q =a.P. We translate the process by including the command 


[ax,y] (st = sq)&(loc = locy) — (st' = sp); 


for each location / and each module y that may perform an input on channel 
a. We then append the translation of P. 


Case 4: Q = \/.P. We translate the process by including the command 
[tick] (st = sq) — (st’ = sp); 
and appending the translation of P. 
Case 5: Q = rep,;.P. We translate the process by including the command 
[repy.2] (st =8q) —> (st = sp)&(loc’ = locy); 


for each module y that may output on channel rep. We then append the 
translation of P. 
Case 6: Q = Teéps.P. We translate the process by including the commands 
[| (st = sg) &(i,s > 0) — (i, =i, — 1)&(st’ = stg, ); 
[repz.y] (st = stgo) —+ (st' = sp); 


for each inactive module y of species s, and then appending the translation 
of P. 


Case 7: Q = a1.P, + a2.P2. We translate the process by computing the 


translations [a1.P;] and [a2.P2] and replacing all commands of the form 


[act] (st = sta,.p, )&guard — updates; 
by 


[act] (st = stg)&guard —> updates 
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and similarly for the commands of [a2.P.]. 


Case 8: Q=p,: Pi+...+ pn: Py. We translate the process by appending 
[Pi],.--, [Pr] to the command: 


I] (st = sq) —> pi: (st’ = sp,) +...+-pn: (st’ = sp, ); 


Case 9: Q = cond (e; > Pi,...,€n & Py). We translate the process by 
constructing [Pi],..., [Pn] and replacing each command of the form: 


[act] (st = sp, )& guard — updates; 
by the command 
[act] (st = sQ)&[e1 Gloc] &...&e:—1 Gloc] &e; @loc]& guard —> updates; 
where [eQ/] is the translation of the PALPS expression e@f into the PRISM 
language. Note that although all the transitions are simultaneously enabled 
from state sg, the use of the condition assigned to each of the conditional 


alternatives implies that at most one P; can be selected and this is such that 
i is the smallest 7 with e; evaluating to true. 


Case 10: C, C te! P. We translate the process by computing [P] and 
replacing each command in [P] of the form 

[| (st = sp)&guard — updates; 
by 


[| (st = sc)& guard — updates; 


Case 11: Q =0. We translate the process as 


l] (st = so)&(loc = 1) —> (s} = 8 1)&(st! = Sinactive)} 


[tick] (st = Sinactive) —— (st = Sinactive) 


The translation above is mostly self-explanatory. It is worth pointing out, 
however, that, in Case 6, the translation of the activation on a new individ- 
ual is conducted in two steps: in the first step, it is confirmed that there 
is still an inactive individual to be activated, (is > 0) before reserving this 
individual via the assignment i, = i, — 1 and, in the second step, module x 
synchronizes with one such individual y via action rep, y. 
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4.4 Correctness of the Translation 


We now turn to consider the correctness of the proposed translation. This is 
demonstrated via the following two theorems. In what follows, given a PRISM 
model M, we write M —+ M’ if M contains an action [a] guard -> updates; 
where guard is satisfied in model M and execution of the updates gives rise 

to model M’. 


Theorem 1. For any configuration (£, Sys), where EF is compatible with 
Sys, whenever (E, Sys) “+ (E’, Sys’) then either [(E, Sys)] —> [(E’, Sys’)] 
or [(E, Sys)] —— [E’, Sys’)]. 


Theorem 2. For any configuration (£, Sys), where EF is compatible with 
Sys, whenever [(E, Sys)] —> M then (E, Sys) > (E', Sys’) and one of 
the following holds: 


1. either M = [(E£’, Sys’)], or 


2, M — M’, M’ = [(E£’, Sys’)] and for all M —>+ Mj, we have 
(E", Sys’) aes (E"", Sys”) and M,; —> Mp2 where Mp = [(E”, Sys")]. 


Theorem 1 establishes that each transition of (FE, Sys) can be mimicked 
in one or two steps by its translation module. Theorem 2 considers the other 
direction of the correctness: Given a transition of a PRISM module there are 
two possibilities. First, it is possible that the transition of the module corre- 
sponds to exactly one transition of the associated PALPS process. Second, it 
is possible that the transition of the module has resulted in an intermediate 
step of establishing the transition of a PALPS process. In this case, according 
to the theorem, each transition M —> M, of this intermediate state M can 


also be executed on the side of the PALPS process, (E’, Sys’) asa (E"", Sys"), 
and the pending move of the first transition can still be executed yielding a 
state Mz which is in fact the translation of (E”, Sys”). 


Sketch of the proof of Theorem 1: The proof consists of a case 
analysis of all possible ways in which the transition (EF, Sys) “> (E’, Sys’) 
can be produced. Three cases exist: 


e If the transition involves a single process participant P:(s, ?), then we 
may verify by induction on the structure of P that any action of P can 
also be performed by its translation and the resulting PRISM model is 
the translation of (E’, Sys’). In all cases this can be established in a 
single move of module P. 
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e If instead the transition has arisen via the communication of two com- 
ponents of Sys then it is possible to establish that the two modules 
corresponding to the two components share the action in question and 
can thus execute the synchronization. This may take one or two ac- 
tions at the PRISM level depending on whether process activation is 
involved. 


e Finally, if a = \/, then it must be that all components of Sys enable 
the transition \/. We may then observe that the PRISM translations 
of the components enable the tick action and thus the transition can 
be performed in exactly one move. In any case, the resulting PRISM 
model is the translation of (E’, Sys’). 


Sketch of the proof of Theorem 2: The proof consists of a case 
analysis of all possible ways in which the transition [(£,Sys)] —> M can 
be produced. It follows along similar lines with the proof of Theorem 1, with 
the exception of the first step implementing the activation of an inactive 
module. The important point to note here is that the intermediate state MW 
captures correctly environment E’ in the transition (EL, Sys) —> (E’, Sys’). 
Thus, any further transition from this intermediate state M is in fact also 
enabled at the PALPS level. The resulting states are equivalent, pending the 
second step of the activation of the inactive individual. 


4.5 Verification in PRISM 


In this section we briefly describe the types of analysis that can be performed 
on PALPS models via the PRISM model checker. For details, we refer the 
reader to [8]. 


Model Checking To begin with, PALPS models may be model checked 
in PRISM against properties specified in the PCTL logic. The syntax of the 
PCTL logic is given by the following grammar where ® and @ range over 
PCTL state and path formulas, respectively, p € [0,1] and k € N. 


® true | e | -® | ®A® | Propld] 
bd := XO | GUS | O,UG 


In the syntax above, we distinguish between state formulas ® and path 
formulas ¢, which are evaluated over states and paths, respectively. A 
state formula is built over logical expressions e and the construct Py.p[¢]. 
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Intuitively, a configuration s satisfies property P..»|¢] if for any possible 
execution beginning at the configuration, the probability of taking a path 
that satisfies the path formula ¢ satisfies the condition bX p. 

Path formulas include the X (next), U* (bounded until) and U (until) 
operators, which are standard in temporal logics. Intuitively, X® is satisfied 
in a path if the next state satisfies path formula ®. Formula ®; U'®, is 
satisfied in a path if ®; is satisfied continuously on the path until ®2] becomes 
true within & time units steps. Finally, formula ©,;U®, is satisfied if ®2 is 
satisfied at some point in the future and ®, holds up until then. 

As an example, consider a population s in danger of extinction. A 
property that one might want to check for such a population is that the 
probability of extinction of the population in the next ten years is less than 
a certain threshold pe. This can be expressed in PCTL by the property 
Pcp, [trueu? pep 5¢8@ = 0]. Alternatively, one might express that a 
certain central location will be re-inhabited with at least some probability 
pr by sQ@l = 0 > Psy, [trueU(s@é > 0)]. 

It is also possible to study the relation within a model between the size 
of the initial population and the probability of extinction of the population, 
by checking properties of the form s@é > m — Psp, |trueU(s@é = 0)] or to 
explore the dynamics between two (or more) competing populations s and 
s’. As an example, expressing that, within the next 20 years with some high 
probability, members of the population s will outnumber the members of 
population s’: Psp[trueU?9(S>pepog 8’ QL < Vvetroe S@2)]- 

PRISM also enables to take a more quantitative approach for model 
checking PCTL properties: it supports the verification of the constructs 
Prin—?|¢) and Pmar—2?|¢| via which the minimum and maximum probabili- 
ties of satisfying ¢ are computed. 


Steady-state behavior. PRISM also supports reasoning about the steady- 
state behavior of a model, that is, the behavior in the long-run or when 
an equilibrium is reached. Steady-state properties are only available for 
discrete-time and continuous-time Markov chains. These properties are ex- 
pressed by Sjoung[prop|. Such a property is true in a state s of a discrete- 
time or a continuous-time Markov chain if, starting from s, the steady-state 
(long-run) probability of being in a state which satisfies the (boolean-valued) 
property prop, meets the bound bound. For example, the steady-state prop- 
erty Sbouna|s@2 = 4] expresses that the long-run probability that there will 
be 4 individuals of species s at location 2 meets the bound. 
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Rewards. PRISM models can also be augmented with information about 
rewards: It is possible to assign a reward (a positive real number) to any 
command or state of a PRISM model. Every time a command is executed 
or a state is visited, the rewards associate with the command or state is 
accumulated. It is then possible to reason about reward-based properties for 
discrete-time Markov chains, by extending the logic PCTL with the following 
additional operators [18]: 


Roar [C"] | Roor[IS"] | Roar[F®] | Roar |S] 


where ie {<,<,>,>},r € Rso,k € N and ® is a PCTL formula. The R 
operator defines properties about the expected value of rewards. The formula 
Rear |W], where = denotes one of the four possible operators in the grammar 
above, is satisfied in a state s if, from s, the expected value of reward w 
meets the bound = r. Operator C<* refers to the reward accumulated 
over k time steps; I<* the state reward at time instant k; F®, the reward 
accumulated before a state satisfying ® is reached; and S, the long-run rate 
of reward accumulation. Properties of the form R=?[~] means “what is the 
expected reward for operator ~?”. 


5 A Case Study in PRISM 


In this section, we apply our methodology for model checking PALPS sys- 
tems using the PRISM tool and we explore the size of models on which 
model checking can be performed. As an example we consider a number of 
variations of the system in Example 1, described in Section 2. 

In our model we will assume a lattice of locations of size n x n (for 
n = 2,3,4). To begin with, we considered a lattice of 2 x 2 with reflecting 
boundary conditions. This implies that, in case of migration, an individual 
may move to one of two other locations. Thus, the PALPS definition of an 
individual takes the following form: 


1 
pe »~ 5° gol.,/.cond (sQmyloc = 1 > Pi; true> /.P) 
LENb(myloc) 
def 

Pi FS p:rep./.P, ® (1— p) : rep rep./.P, 

The PRISM encoding of the system follows the translation methods pre- 
sented in Section 4.3. It contains four global variables, s1, s2, s3, s4, which 
represent the number of individuals at each location. It also includes the 
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constant p which characterizes the probability that an individual produces 
one or two children. The PRISM model contains a module for each active 
individual and a module for each inactive individual. All the tests were 
performed on a G46VW Asus laptop with an Intel i5 2.50 GHz processor 
and 8 GB of RAM. We ran the tests under Linux Ubuntu 13.04 (Kernel 
3.8.0_17), using PRISM 4.0.3 (with the MTBDD engine) and Java 7. 

As a first experiment we explored the size of a system that can be 
handled for model checking in PRISM by setting a limit on the maximum 
number of individuals that can be created. It turns out that the size of a 
model with seven or more individuals exceeds the amount of RAM memory 
available. Furthermore, we created a variation of the model with one preda- 
tor in addition to the individuals of type P above. In this case it turns out 
that the size of the model with 1 predator and 5 preys, exceeds the size 
that can be treated in PRISM on personal computers which, according to 
the PRISM manual, is 10’ — 108 states. 


Case study Number of | Number of | Construction RAM 
size States Transitions | time (sec.) (GB) 
3 PALPS individuals 130397 404734 8 0.5 
2 PALPS individuals* 148360 250964 13 0.8 
4 PALPS individuals 210736 7312132 101 1.9 
5 PALPS individuals 18647917 85022112 601 2.8 
6 PALPS individuals | 87332123 | 116457321 5094 4.5 
4 PALPS individuals* | 113450511 | 411434511 3521 5.1 
7 PALPS individuals out of RAM 


Table 5: Performance of building probabilistic models in PRISM. 


In Table 5 we summarize the results we obtained in our experiments. In 
the models, we fixed p = 0.4. Results marked by * refer to models extended 
with a predator. 

We believe that the reason for the state-space explosion is the large 
number of possible interleaving orders among all the commands of all the 
modules, which, in the case of model checking have to be computed in their 
totality. In next section we discuss how to overcome this problem in the 
future. 

As a second experiment, we proceeded to determine the limits for sim- 
ulating PALPS models. We constructed PRISM models with various numbers 
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of modules of active and inactive individuals and we run them on PRISM. In 
Table 6, we summarize the results. It turns out that for models with more 
than 5000 individuals simulation requires at least 12 hours (which was the 
time limit we set for our simulations). 


Individuals | File Size (MB) | RAM (GB) | Simulation Time (s) 

10 0.1 0.18 1 

100 0.4 0.3 8 

500 2.0 0.5 45 

1000 4.2 1.0 300 

1500 6.2 0.7 454 

2000 8.2 0.9 820 

5000 20.1 2.0 > 12 hours 

10000 44.1 3.4 > 12 hours 


Table 6: Performance of simulating probabilistic systems in PRISM. 


Subsequently, we looked into the restriction imposed by our assumption 
of bounded replication. This restriction may lead to deadlocks when an 
active individual attempts to reproduce but no inactive module is available 
for synchronization. To explore this, we used the simulation environment 
of PRISM and searched for deadlocks by repeating 100 simulations of the 
model of maximum path length 1000 time steps. Although, this procedure 
is not complete, we may consider it sufficient as it looks into a fairly large 
number of life cycles of the population. In Table 7 we summarize the results 
obtained. 

As a further experiment with PRISM, we considered a simplification 
of Example 1 in which the model consists of exactly 3 individuals and in 
which reproduction does not produce new processes but only updates the 
number of individuals at the reproduction location. Therefore, there will 
always be 3 active modules, but the population size at each location will be 
updated every time an individual performs the reproduction action. The 
PRISM model of this example contains 4 modules, one for each of the three 
individual, and one module for modeling reproduction. As an example, we 
present some excerpts of the module to define one of the individuals, namely 
Pi. First, the individual has to choose to move to a new location. This 
is modeled as follows where the variables st1 and loci represent the local 
state of the module and the location of the individual, respectively. 
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Active individuals | Inactive individuals | Deadlock (2, y, z) 
3 18 (No, Yes, Yes) 
4 24 (No, Yes, Yes) 
5 30 (No, Yes, Yes) 
6 36 (No, Yes, Yes) 
7 42 (No,No, Yes) 
8 48 (No,No, Yes) 
9 56 (No,No, Yes) 
10 60 (No,No,No) 


Table 7: Occurrence of deadlock in various instances of the model. The 
values (x,y,z) refer to the presence of deadlock in the case of 4 locations 
(x), 9 locations (y), and 16 locations (z). 


sti : [0..40] init 0; 
// For location 1 
(] (st1=0)&(loc1=1) -> 
0.5: (st1’=17)& (si’=s1-1)&(s2’=s2+1)&(loc1’=2) 
+ 0.5: (st1’?=27)&(s1’=s1-1)&(s3’=s3+1)&(loc1’=3) ; 


Subsequently, the process checks if it has absolute control of the location. If 
it does, then it reproduces; otherwise, it synchronizes with the clock action. 


// For location 1 
[] (sti = 17) & (loci=1) & (st = 1) -> (st1’ = 19); 
[] (sti = 17) & (loci=1) & (sit != 1) -> (st1’ = 15); 


In the case that reproduction is selected, the individual chooses to reproduce 
one or two times. 


// For location 1 

[] sti = 19 -> pp: (sti? = 12) + (1-p) : (sti’ = 11); 
[rep1] (st1 =11) -> (st1’ = 12); 

[rep1] (st1 = 12) -> (st1’ = 15); 


Whether the individual reproduces or not, it synchronizes with the other 
modules with a clock action and, finally, returns to the initial state. 


[tick] (sti = 15) -> (sti’ = 0); 


158 A. Philippou, M. Toro, M. Antonaki 


The reproduction module is defined as follows. Local variable st4 rep- 
resents the local state of the process. At state 0, the module can synchronize 
with the other modules via the clock action or the module can synchronize 
with any individual at any location with the respective reproduction action. 
After synchronizing with the reproduction action, the module goes to a state 
in which it updates the number of individuals at the respective location by 
one and then returns to the initial state. 


st4 : [0..34] init 0; 

[tick] (st4 = 0) -> (st4’=0); 

// For individual 1 

[rep1] (st4=0) -> (st4’=1); 

[] (st4=1) & (loct=1) -> (st4’=0) & (si’=s1+1); 
[] (st4=1) & (loct=2) -> (st4’=0) & (s2’=s2t+1); 
[] (st4=1) & (loc1=3) -> (st4’=0) & (s3’=s3+1); 
[] (st4=1) & (loct=4) -> (st4’=0) & (s4’=s4+1); 


As a first analysis of this simplified example, we studied the effect of 
the value of the probability p on the size of the population by performing a 
reachability analysis of the model. Specifically, Figure 5 records the proba- 
bility to reach a state with a certain total number of individuals for various 
values of the constant p. 


Probability to reach a state with n individuals 


Probability 


Probability to 
reproduce (p) 


is Individuals(r) 7 
Figure 5: Probability that the population reaches a certain size, for various 
values of p. 


Using PRISM it is also possible to study the steady-state behavior of 
the model. For instance, one might wonder about the effect of the initial 
placement of individuals on the grid for the long-run behavior of the sys- 
tem. To this effect, we have enunciated and checked the following property 
S =? [s, = r], where q € [1,4] is the identifier of the location and r € [0, 6] 
the number of individuals. This property allows us to evaluate the proba- 
bility to have r individuals at location qg in the long run. For example, for 
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q = 1,r = 2, this property means “what is the probability that, in the 
long run, we reach a state in which the number of individuals at location 1 
is equal to 2?”. Interestingly, our experiments have shown that the initial 
placement of individuals has no effect on the long-run behavior of the model: 
in the long run all locations have the same probability of being populated 
by a population of r individuals for any value of r. Intuitively, this is so 
because as there are no predators in the system, in the long run the pop- 
ulation densities of the locations balance out. In Table 8 we present these 
probabilities for r € [0, 6]. 


| Number of individuals | 0 | 1 2 3 4 a. | 6x] 
| Probability 0.0 | 0.191 | 0.418 | 0.261 | 0.109 | 0.021 | 0.0 | 


Table 8: Steady-state probability for reaching a state with a certain number 
of individuals at a location 


Finally, in Table 9, we present some results that we obtained by as- 
signing rewards to certain actions of the model. In particular we assigned 
rewards to (1) the clock action and (2) the reproduction actions. We then 
used PRISM with the property R =?[C’S*] to compute the average reward 
accumulated within the first k execution steps. In the case that rewards are 
associated with the clock action, this property calculates the average num- 
ber of clock steps taken within the first k execution steps. Similarly, in the 
case that rewards are associated with the reproduction action, this property 
calculates the average number of reproductions within the first k execution 
steps. In the latter case, we observed that the results were the same for each 
individual. That is, each individual engages in the same average number 
of reproductions within the first & execution steps. Intuitively, this is so 
because there are no predators in the system and, following the first move 
action and thereafter, each individual is equally likely to reproduce. 

As an example of how this is implemented in PRISM, to count the 
average size of the offspring of individual P1, we add a reward of 1 each time 
individual P1 synchronizes with the reproduction module by a reproduction 
action. This reward structure is defined as follows: 
rewards "repPi" 

[repi] true : 1; 


endrewards 


Thus, as reported in Table 9, the expected number of tick-synchronizations 
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within the first 100 steps of the model is 10.558, whereas the expected num- 
ber of reproductions an individual engages in during the first 100 steps of 
the model is 2.117. 


Simulation steps | Expected time units | Expected size of offspring 
100 10.558 2.117 
1000 139.049 2.133 


Table 9: Expected number of clock steps and the expected number of repro- 
ductions of an individual 


6 Concluding Remarks and Future Work 


This paper reports on work towards the development of a process-calculus 
framework for the spatially-explicit and individual-based modeling of eco- 
logical systems which follows a discrete-time, discrete-space approach. The 
main distinguishing feature of PALPS is that it associates locations with at- 
tributes which allow the enunciation of location-dependent behavior. For 
example, “an individual will migrate from its current location only if this 
location has exceeded its capacity”. This is achieved by enriching the syn- 
tax with a set of expressions and extending the notion of a process by the 
notion of a configuration. A configuration is an entity that combines the 
process/code of a system with an environment which records information 
regarding the state of the system, based on which expressions are evaluated. 

Furthermore, we have described a methodology for translating PALPS 
models into the PRISM language. The proposed encoding can be employed 
for both simulating and model checking PALPS systems using the PRISM tool. 
In this paper we experimented with the model-checking capability and we 
demonstrated types of model-checking analysis that can be performed on 
PALPS models via simple examples. This exercise has confirmed that model 
checking of population systems defined in PALPS is not viable, at least at 
this stage. This is because, the large number of modules and the high degree 
of nondeterminism that occurs within PRISM encodings of PALPS modules 
results in a state-space explosion that becomes apparent even in systems 
with a very small number of individuals. 

As far as simulation is concerned, we have observed that PRISM allows 
the simulation of realistic PALPS models. Nonetheless, the individual-based 
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approach implemented by PALPS and its translation into PRISM impose lim- 
its on the size of systems that can be simulated that are not as prominent 
in stochastic models: As we discussed in Section 1, it is indeed the case that 
continuous time frameworks can benefit from Gillespie’s stochastic simula- 
tion methods, which speed up discrete stochastic simulation by leaping over 
many reaction events in a way that approximates exact stochastic simula- 
tion. 


These conclusions lead to the challenge of defining new approaches for 
analysis of MDPs that arise from the modeling of population systems. A 
number of directions may be employed for tackling this challenge. The 
first method is related with the idea of sequentializing execution of mod- 
els by imposing an ordering between actions and processes. As we have 
already discussed, this idea has stemmed from ecologists concerned with 
the construction of discrete-time ecological models: According to [5], in 
discrete-time, individual-based models, ordering of processes and actions 
can be applied on models and different policies on these orderings can be 
explored and studied separately. 


Applying such orderings on PALPS would greatly simplify the state- 
space of models as multiple interleaving orders that are causing the models 
to explode will be pruned away and nondeterminism will be significantly 
reduced. We also expect that sequentializing the execution of different types 
of actions will enable the grouping of sequences of actions into single macro- 
actions thus yielding smaller models while preserving the initial behaviors. 
For example, given a sequence of movement actions by a set of individuals we 
can simply retain the first and last state while pruning away the intermediate 
states of the separate individuals moves. We are currently exploring this 
approach by enriching PALPS models with policies that determine externally 
the orderings of different actions within a model. To this end, we are also 
interested in developing the notion of (partial) confluence [39] for PALPS and 
applying it for reducing the state space of PALPS systems. 


Other possible approaches we plan to investigate for reducing the state- 
space of PALPS models include enhancing the semantics of PALPS to enable a 
more succinct presentation of systems especially in terms of the multiplicity 
of individuals. Such a semantics would have as its basic entity the process 
P:(s,¢,n), which represents n individuals of species s at location ¢, and 
it would ease both the specification as well as the verification of system 
models. In addition, we are interested to explore the possibility of giving 
symbolic semantics on PALPS systems by applying a symbolic treatment of 
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environments. Furthermore, an interesting alternative would be to extend 
the work of [30] in the spatial domain for defining mean-field semantics for 
PALPS. 


As far as the PRISM tool is concerned, we believe that it would be 
worth to consider alternative approaches for producing PRISM input from 
PALPS models. We recall that there are two possibilities for achieving this: 
using the PRISM language or constructing the transition matrix of a Markov 
decision process. Adopting the second approach could prove beneficial for 
a number of reasons: (1) In PALPS, we use CCS communication instead of 
the CSP communication employed in PRISM language. This is necessary for 
modeling phenomena such reproduction and preying where the communi- 
cation is binary and retaining the asymmetry between the roles of input 
and output is important for implementing the task (e.g. we must distin- 
guish between the actions of a prey and a predator). (2) Direct encoding 
of PALPS into a transition matrix would avoid the bounded replication and 
subsequent modeling of inactive individuals in a PRISM model, deemed nec- 
essary due to PRISM not allowing the dynamic creation of modules. Note 
that the translation of recursion for the probabilistic pi-calculus into PRISM 
proposed in [34] is not sufficient for our purposes because the language con- 
sidered for the translation concerns a restricted kind of recursion and no 
replication. (3) Finally, another source of complexity in the translation con- 
cerns the modeling of reproduction, where an intermediate step is required 
to model the process due to the necessity of updating a global variable in a 
synchronization command, something that is not allowed within the PRISM 
language. These observations suggest that representing PALPS models as a 
Markov decision process, instead of encoding them into the PRISM language, 
could result in simpler models. 


Another possible direction for future work which aims to enhance the 
expressiveness of PALPS for capturing more complex systems includes the 
adoption of continuous time as well as the use of dynamic attributes to 
allow exploring the system while, for example, patch quality degrades, tem- 
peratures increase, etc. This could be made possible via the combination 
of continuous and discrete time behaviors similarly to the behavioral hybrid 
process calculus of [28]. Finally, we are currently involved in the application 
of our methodology to case studies from the local habitat, the intention 
being to explore properties such as species extinction and persistence. 
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